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Abstract 

Oscillatory chemical reactions often serve as a timing clock of cellular processes in 
living cells. The temporal dynamics of protein concentration levels is thus of great 
interest in biology. Here we propose a theoretical framework to analyze the frequency, 
phase and amplitude of oscillatory protein concentrations in gene regulatory networks 
with negative cyclic feedback. We first formulate the analysis framework of oscillation 
profiles based on multivariable harmonic balance. With this framework, the frequency, 
phase and amplitude are obtained analytically in terms of kinetic constants of the reac- 
tions despite the nonlinearity of the dynamics. These results are demonstrated with the 
Pentilator and Hes7 self-repression network, and it is shown that the developed analy- 
sis method indeed predicts the profiles of the oscillations. A distinctive feature of the 
presented result is that the waveform of oscillations is analytically obtained for a broad 
class of biochemical systems. Thus, it is easy to see how the waveform is determined 
from the system's parameters and structures. We present general biological insights 
that are applicable for any gene regulatory networks with negative cyclic feedback. 

Keywords: Gene regulatory network; Periodic Oscillations; Harmonic Balance; Systems bi- 
ology 



1 Introduction 
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In living cells, oscillatory chemical reactions serve as a timing clock of important cellular 
processes. The tempor al dynamics of oscillatory gene expression has thus been actively 
studied in recent years (jGoldbeter and Berridgel 11997 ; IWinfred . l200lh . It is known that 
the frequency, phase and amplitude of the oscillations are diverse, ranging from minutes to 
hours and from phase synchronization to asynchronous oscillations. However, the relation 
between the dynamical properties of biochemical system and the resulting temporal pattern 
is not thoroughly understood. Here we propose a theoretical framework for quantitatively 
studying the frequency, phase and amplitude of oscillatory chemical concentrations that 
arise from biochemical networks with negative cyclic feedback. 

The negative cyclic feedback motif shown in Fig. [TJ where each gene activates or represses 
another gene expression in a cyclic way, ha s been considere d as a core ci r cuit module t o 
produce periodic oscillations for a long time ( Goodwin! ( 19651 ); [Tyson ( 1975F ); Thronl ( 199l[ ); 
Hori et al.l (|201ll) for examp le). This conjecture was re cently corroborated with a synthetic 
biological circuit reported in lElowitz and Leibler in which oscillatory gene expression 

was indeed observed in a gene regulatory network consisting of three repressors. Moreover, 
the negative cyclic feedback was also found in many existing bi ochemical networks that 



exhibit periodic oscillations such as the somitogenesis oscillator (jHirata et al.l . 120041) and 
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p53 networks (jLahav et al.1 . 120041) . Thus, the study of negative cyclic feedback circuits can 



potentially unravel the essential dynamical properties of oscillatory biochemical reactions. 
The frequency profile of oscillatory gene expression in cyclic biochemi cal networks was 



studi ed by means of numerical simulations for the Repress i lator motif (lEl-Samad et al 
20051 ) and for the self-negative fee dback motif dMonkl . 120031: IWang et all 12004) . A more 



theoretical appr oach was t aken by IRapp ( 1976 ) to investigate the fre quency of Go odwin 
type oscillators ( Goodwin . 19651 ) based on harmonic balance analysis (Khalil, 2001). This 
approach is, however, not directly applicable to the biochemical oscillators of our interest 
due to the multiple nonlinearities, and the phase profile of oscillations is still open to be 
analyzed. Hence, a more general theoretical framework is desirable to systematically explore 
the profiles of oscillations in biochemical networks. 

The goal of this paper is twofold: (i) to develop a systematic method to study the profiles 
of oscillations in gene regulatory networks with negative cyclic feedback, and (ii) to gain 
novel biological insight on the relation between the oscillation profiles and the reaction 
kinetics. In particular, we here pursue universal insights that are applicable to a broad class 
of biochemical networks rather than the numerical study of particular biological examples. 
To this end, we adopt a simplified model that is suitable for capturing the essence, and we 
analytically obtain the oscillation profiles. Spe cifically , we use describing fu nction approach 
to approximating nonlinearities of the system (|Khalill . l200ll : Ilwasakill2008l) . and formulate 
multivariable harmonic balance equations for the analysis of biochemical oscillators. Using 
distinctive features of the dynamics of biochemical networks, we derive analytic estimates 
of frequency, phase and amplitude in terms of kinetic constants of the reactions despite the 
nonlinearity of the dynamical model. 

A distinctive feature of our result is that the approximate waveform of oscillations is an- 
alytically obtained for a broad class of cyclic biochemical systems. Thus, it is easy to see 
how the waveform is determined from the system's parameters and structures. This feature 



is dem onstrated through the analysis of two biolog i cal ex amples, the Pentilator (jTsai et al 



20081 ) and the somitogenesis clock ([Hirata et all 120041 ) . The Pentilator is a conceptual 
biochemical oscillator with five repressors connected in a cyclic way. We here analyze the 
Pentilator to show that the developed result is useful even for large-scale gene regulatory 
networks. _In the so mitogenesis example, we investigate the Hes7 self-repression network 
20041 ) . It was experime ntally observed that Hes7 gene exhibits oscillatory 



(Hirata et al 



expression with almost two- hour cycle ( Hirata et al. . 20041 ). We demonstrate that the de- 



veloped analysis method indeed predicts the two-hour period, and show that transcription 
and translation time delays play an important role in maintaining the oscillation period. 
Finally, general biological insights that are applicable to any biochemical networks with 
negative cyclic feedback are presented. 

This paper is organized as follows. In Section[21 the dynamical model of the gene regulatory 
network is introduced. Then, we formulate multivariable harmonic balance for biochemical 
oscillators in Section |3l Based on this formulation, the frequency, phase and amplitude of 
oscillations are analytically obtained in Section 21 Section [5] is devoted to illustrate these 
results with the biological examples. In Section [6j we provide general biological insights. 
Finally, concluding remarks are given in Section [7] 

The basic idea of the ana lysis was previously demonstrated in the authors' conference 
paper ( Hori and Hara , [201ll ) with the omission of rigorous proofs. In this paper, we include 
complete proofs of the theorems. Moreover, the analysis of oscillation amplitude (Section 
I4.3[) and the illustrative examples of the two biochemical oscillators (Section[5]) are presented 
as original work. The biological insights in Section |6] are also greatly extended. 
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Figure 1: Gene regulatory networks with negative cyclic feedback. The symbols — ¥ and H 
represent activation and repression o f transcription, respe c tively . (Left) activator-repressor 
motif, (Center) Repressilator motif (jElowitz and Leiblerl . 120001) . (Right) generic negative 
cyclic motif considered in this paper. 



2 Model Description and Existence of Oscillations 

2.1 Dynamical Model of gene regulatory networks with negative 
cyclic feedback 

In cells, the production of proteins consists of two processes, transcription and translation. 
In transcription, information on a gene is decoded, and turns into the copies of messenger 
RNA (mRNA) . The mRNA molecules are then translated to corresponding protein molecules 
(see Fig. [1]). These proteins then activate or repress the transcription of other genes, and 
form gene regulatory network. 

In this paper, we consider the gene regulatory network illustrated in Fig. [T] Here each 
protein activates or represses another transcription in a cyclic way. The d ynamics of mRNA 
and protein concentrations in su ch networks can be modeled as follows ( Chen and AiharaL 
2002t lElowitz and Leiblerl . l2000h . 



hit) 

Pitt) 



c i r l {t - r n ) - hpi(t), 



(1) 



for i = 1,2, • • • , N, where r, € M + (:= {x 6 R | x > 0}) and pi € R+ denote the con- 
centrations of the i-th mRNA and its corresponding protein synthesized by the i-th gene, 
respectively. For the sake of notational simplicity, we regard the subscript as N through- 
out this paper. The positive constants r ri and r Pi describe transcription and translation 
time delays resulting from unmodeled intermediate process, respectively. The kinetic con- 
stants a,-,6,-,Cj and /3j represent the followings: and bi denote the degradation rates of 
the i-th mRNA and protein, respectively; Cj and /3j denote the translation and transcription 
rates, respectively. The nonlinear function /j(-) : R + — > R + stands for the effect of cither 
activation or repression of the transcription, and it is defined by 
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/i(p*-i):= 



1+PT-i 
Pi-i 

1+P?-! 



— (= : F R (pi _i ) ) (repression) 



(=: FAipi-x)) (activation), 



(2) 



with a Hill coefficient v. 

We see from ([1]) that the transfer function from u% := fi(pi-\) to pi, which describes 
the dynamics of protein production in each gene, becomes a second order system with time 
delay. Defining these second order systems as hi(s) (i = 1,2, ••• >N), the overall system 
dl]) can be described by Fig. [2] (Left), where a transfer matrix H(s) and a static vector 
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Figure 2: (Left) The block diagram of the gene regulatory networks with negative cyclic 
feedback, where p := \pi(t — T pl ),p2(t — T P2 ),--- ,pN(t — T PN )] T . (Right) linear system (s) 
in (|16|) . The static nonlinearity / is replaced with the corresponding describing function. 



nonlinearity function f are defined as 

H(s) := diag(fci(«), h 2 (s), • • • , h N (s)), (3) 
/:=[J^/i(0,-^/a(0»--- ,R 2 N fN(-)] T - (4) 
The transfer function hi(s) and the constants Ri are 

e -s(r ri +T p .) ^ 1 

(T Qi s + ±){T bi s + 1) a, 6i 

fl .. = vg|. ( i = l,2,---,iV). (6) 



It should be emphasized that i?(s) describes the dynamics of each gene's protein production, 
and / describes the interaction between genes. 

The dynamical behavior of the system (fT|) can be qualitatively classified by 



TT r r 1+1 (/<(0 = *U(0) ^ 



Specifically, the protein concentrations asymptotically converge to one of equilibr ia when 
<5 > 0, while they exhibit oscillatory solutions as well as convergence when d < ( Encisol 



l2007tlSmithlll987h . Therefore, the gene regulatory networks with S < are of interest in 
this paper. 

Assumption 1. For given /,(■) (i = 1, 2, • • • , N), 8 < 0. 

This assumption implies that a given gene regulatory network has an odd number of repres- 
sive interactions (dfi/dp < 0) between genes. Thus, the loop gain of the overall system is 
negative. 

In order to capture the essential dynamical properties in an analytic way, we hereafter 
consider a simplified model of (JJ). It is assumed that the kinetic parameters a,, Cj and ft 
in (JTJ) are homogeneous between genes throughout this paper. 

Assumption 2. a\ = a 2 = • • • = aAr(=: a), 6i = 62 = • • • = &jv( =: c i = c 2 = ' ' ' — 
c N (=: c) and ft = ft, = •• • = ftv(:= /?) in ©. 

2.2 Existence of oscillations 



In iTakada et~ai1 (|20ld) . the authors presented existence conditions of oscillations for the 



gene regulatory system ([T]). The analysis then revealed that five dimensionless quantities 
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Table 1: Physical meanings of the constants 



N The number of genes in gene regulatory network 
Q Discrepancy of mRNA and protein degradation 
time 

R Ratio of degradation and production rates, 
which accounts for equilibrium concentrations 

f Mean time delay in transcription and transla- 
tion normalized by the mean degradation rate 

v Hill coefficient, which quantifies the degree of 
cooperative binding 



essentially determine the existence of oscillations, namely (AT, Q, R, f , v), where 

Q-=t^, # : =-r> T '-=Tf- 

T A ab T A 

with 



T G c/3 _ t T a + T b 

R: =zt:> r: =^ t a ■= — o — ( 8 ) 



T a .= -, T b .= -, T.= . (9) 

Physical meanings of these constants are summarized in Table [TJ 

Remark 1. Since Q is the ratio of geometric and arithmetic means, < Q < 1 holds. In 
particular, Q — > 1 implies that the degradation time constants of mRNA and proteins, T a 
and T(,, tend to be a same value. The parameter f describes the ratio of mean transcription 
and translation delay to the mean degradation rate. 

In what follows, we analyze the frequency, phase and amplitude of oscillatory protein 
concentrations, and reveal how these parameters in Table Q] affect the oscillation profiles. 
Throughout this paper, it is a s sumed that the system (fTJ) satisfies the existence conditions 
of oscillations in Takada et al. ( 2010l ). and exhibits oscillatory protein concentrations. 



3 Multivariable Harmonic Balance Analysis for Biochem- 
ical Oscillators 

In this section, we provide a theoretic framework of oscillatio n profile analysis for the gen e 
regulatory network {T]). Using the harmonic balance analysis (|Khalil l2001t llwasakil . 120081 ). 
we first derive a quasi-linear system associated with (TTJ). Then, it is shown that the frequency, 
phase and amplitude profiles are obtained from the equations that the closed loop quasi- 
linear system should satisfy. 



3.1 Multivariable harmonic balance analysis 

The oscillatory waveform of Pi(t) (i = 1, 2, • • • , N) can be written with the infinite sum 
of sinusoidal waves as Pi(t) — J^'kLo aik s ^ n (k^t + tfif.) with some constants w,ipik and 
otik (k = 0, 1, • • • ). Here, we can expect that Pi(t) is approximately written as 

Pi (t) ^Xi+yismizut + ipi) [i = 1,2,- •• ,A r ), (10) 

since the higher order harmonic components are attenuated by the second order low-pass 
filters hi(s) in the cyclic network. It should be noted that Xi and j/i denote the bias and the 
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amplitude of the first order harmonic components of the i-th protein pi (t) , respectively. We 
assume Xi > and yi > without loss of generality. 
Let the describing functions of R 2 fi(-) be defined by 



rn(xi-i,yi-i):= 



2-KXi-! 



fi (xi-i+yi-i sin(t)) dt, 

7T 

fi (xi-i+yi-i sin(t)) sin(t)di. 



(11) 
(12) 



The describing functions rii(xi—i,yi—x) and £j(:Ej_i, yi—i) represent the bias and har monic 
gains of Rffi(-) for the sinusoidal input sin(rot), respectively ( Khalill . 200lh . 

Approximating the nonlinearity /(■) with the describing functions, we see that the input 
to hi(s) is Ui = r)i(xi-i, Jft-i) and it, = &(a;,'_i, J/t-i) instead of Ui = /»(•). Thus, the system 
depicted in Fig. [2] (Left) can be redrawn as shown in Fig. [2] (Right), where 



£o(x, \y\) := cyc(?7i, 772, • • • , r) N ), 
Ki{x,\y\) :=cyc(fi,£ 2 ,--- ,€jv)> 



(13) 



and 



Cyc(Zi,Z 2 ,Z3,--- ,2jv) 



The symbols a; and y are defined as a; := [xi, x%, ■ ■ ■ 
C N with 

tfii := ifi — w 

and \y\ stands for elementwise absolute values, i.e 












• • • Zl 


Z2 








'■■ 




















ZN 



andy := [yie^ 1 , y 2 e ] 



,yNe 



3<?N]T g 



(14) 



tN 



y \ = [2/1, 2/2, 

We can see that the system in Fig. [2] (Right) would satisfy the following closed-loop 
equations, if the waveform of Pi(t) was strictly the sinusoidal wave Xi + yi sm(wt + </?;)• 

(I-H(0)lC {x,\y\))x = (15a) 
(I-H(jm)IC 1 (x,\y\))y = 0. (15b) 

Thus, the solution (vo.x.y) of the equations (fT5|) is expected to capture an approximate 
profile of the oscillations when Pi(t) is sufficiently close to the sinusoidal of the form (fT0|) . 
Consequently, the oscillation profile analysis reduces to the problem of finding 3./V variables 
(za,xi,xt, ■ ■ ■ ,xn, 2/1,2/2! • • • ,Vn, <£2,<P3, • • • ,<Pn) satisfying ((T5J). Note that ipx is defined 
as (fi :— without loss of generality. 



Let x* = [x\,X2,- 



•■N\ 



and y * 



iN 



denote 



the constant vectors that satisfy both of (|T5]) . We define the linear time-invariant systems 
Ho{s) and Hi(s) as 

H.(s) := (I-H^K.y 1 (. = 0,1), (16) 

where K, is the constant matrices defined by K, := K,,(x*,\y*\) (• = 0,1). Although 
there may exist multiple solutions for the equations (|15[) . o rbitally unstable solutions ca n be 
empirically ruled out by checking the stability of Hm(s) (see lKhalil ( 2001 ); Iwasakil (|2008 ) and 
references therein). Specifically, a pair of poles of T-L,{s) is expected to lie on the imaginary 
axis, and the rest in the open left half complex plane, if a solution (w, x*, y*) is orbitally 
stable. Thus, the profiles of stable oscillations can be specified from the solutions of (fTS")) 
and the marginal stability condition. 
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3.2 Oscillation profile analysis for biochemical oscillators 

Using the harmonic balance approach, the oscillation profile analysis of biochemical oscilla- 
tors can be summarized as follows. 



Oscillation profile analysis: Consider the gene regulatory network modeled by ([T]). 
Suppose Pi(t) is approximately written as Pi(t) ~ Xi+yt sm(wt+ipi) (i — 1, 2, • • • , N). 
Then, we find (vcr,ipi,Xi,yi) such that the following two conditions hold. 

(CI) The equations (|T5a|) and (IT5b| hold. 

(C2) The system Ho(s)('Hi(s)) has a pole at s = 0(s = ±jw), and the rest in 
the open left half complex plane. 



In what follows, we explore the solutions of (115]) in an analytic way. We show that the 
equations (|15[) can be rewritten in the form of eigen-equations, and the profiles of oscillations 
can be determined from the eigenvalues/eigenvectors of matrices with a certain structure. 
Let a N x N diagonal transfer matrix U(s) and a scalar transfer function h(s) be defined 

by 

U(s) := diag(e s(T - Tl) , e s(r - T2) , • • • , e s{r - TN] ), (17) 

his) := - t-, r, (18) 

V ; {T a s + l)(T b s + iy 1 ' 



where 



n :=T n +T Pi (< = 1,2,-.. ,N). (19) 



It follows that H(s) — h(s)e~ ST U(s) with r defined in ([9]). Note that H(s) is decomposed 
into the common dynamics h(s)e~ ST and the deviation U(s). Dividing (|15p by /i(0) and 
h{ivj)e^ 3zn and substituting the solution (za,x*,y*), we have 

(0(0)7 - K )x* = 0, (20a) 

(<t)(jm)e jTUT I - UK 1 )y*=0. (20b) 

where 

<f>(s) := l/h(s) and U := U{jm). (21) 

We see that x* and y* can be seen as the eigenvectors of Kq and UKi associated with 
the eigenvalues c/>(0) and 4 , {jw)e^ T , respectively. However, the solution is not obtained by 
simply analyzing the eigenvalue and the associated eigenvector. This is because the matrices 
K a .Ki and U depend on x, \y\ and w as shown in (fT3)) and (|2T|) . thus we need to determine 
both the matrices fC,(x, \y\) and the corresponding eigenvalues/eigenvectors simultaneously. 

Nevertheless, this observation provides us with important intuitions for our subsequent 
analysis. In the following sections, we show that (|20|) can be analytically solved by using 
the structure of the matrices fC,(x, \y\). 

4 Profiles of Oscillations 

In this section, we analytically derive the profiles of the oscillatory protein concentrations 
in terms of the biological parameters in Table [T] 
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4.1 The frequency of oscillations 

We first obtain the frequency w in terms of the parameters shown in Table [TJ The following 
lemma describes the eigenvalue distribution of U{jzu)K\{x, \y\). 

Lemma 1. For any given (m,x,y), the eigenvalues A^ (i — 1,2,- •• ,N) of the matrix 
U(jw)JCi(x, \y\) are given by 



Ai := 



N 



3 J — 7r . (22) 
Proof. For any given (m,x,y), the characteristic equation of U(jvj)lCi(x, \y\) is written 



as 



•s 



N 

N 



n^-i^-iHO. (23) 



Note that the time delays are cancelled out. It follows that Yli=i £( x i-i7 Ui-i) < 0) because 
^i(xi-i, yi-i) is negative/positive when Si is negative/positive, and Assumption 1 holds. 
Thus, ([22]) is derived by solving (gSJ). □ 

We see that the eigenvalues are equiangularly located on a circle with center at the origin. 
Moreover, w does not affect the eigenvalue, and x and y affect only the radial position, but 
not the the angular position of the eigenvalues. Since the oscillation frequency w should 
satisfy (|20|) . possible candidates of w can be characterized as 



vj S 



where arg(-) stands for the argument of a complex number. Note that does not depend 
on x and y. The possible solutions of w are illustrated in Fig. [3l We see from Fig. [3] that 
there are countably infinite solutions. 

Despite the infinite candidates of w, we can show that an orbitally stable solution is unique 
by considering the condition (C2). The following lemma gives a necessary and sufficient 
condition for marginal stability of T-L\{s). 

Lemma 2. The system T-L\(s) defined by (fT6|) has a pair of poles at s = ±juJo and the rest 
in the open left half complex plane, if and only if a pair of the eigenvalues of UKi, which 
we denote by \e and X' e (= Xe), satisfies 4>(joJo) = \t and <p(— jujo) = X' e , and the others lie 
in the domain Cl c + := {7 G C | 4>{s)e ST ^ 7 for Vs G C+}. 

The proof of Lemma 2 is given in Appendix |A] The domain U, c + is the hatched region 
partitioned by <fi(juj)e : > UT in Fig. [3] Lemma 2 implies that the marginal stability of Hi(s) is 
examined from the domain fi^j_ and the eigenvalue distribution of UK\. In particular, the 
gain and phase monotonicity of 4>(s)e ST allows us to show the uniqueness of the frequency 
Vii that satisfies both (CI) and (C2). 

Proposition 1. Suppose there exist (w,x,y) satisfying (CI) and (C2). Then, the fre- 
quency w is uniquely given by the minimum positive solution of aig((j)(j'cu)e :,mT ) = n/N. 

Proof. The distance of the boundary of f2^j_ from the origin, i.e., 



\4>(juW" T \ = V(l + 7>2)(l + T^2 
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• (CI) holds: the candidates of solution 

• (CI) and (C2) hold: the actual solution 









Figure 3: Graphical interpretation of the harmonic balance equation. The red points satisfy 
the bias and the harmonic balance equations, but does not satisfy the marginal stability 
condition. 



monotonically increases for uj > 0. On the other hand, the eigenvalues of UK\ are equally 
distant from the origin as seen in Lemma 1. Thus, the marginal stability condition in Lemma 
2 implies that Ai in (f2"2"|) must lie on the boundary of 51^, when Hi(s) is marginally stable 
(see Fig. [3J. Consequently, w satisfies cf>(jw)e^ T = \\. □ 
From this proposition, we can analytically derive the frequency profile in terms of the 
biological parameters. 

Theorem 1. Consider the gene regulatory networks modeled by (Qp. The frequency zu of 
the oscillatory protein concentrations pi(t)(i = 1,2, ••• ,7V) is expected to be the minimum 
positive solution of 

w = ± (^(^-tnr)+Q 2 -cot(^-tnr)) i. (25) 

Proof. The real and imaginary part of </>(jw) can be written as 

Re^O'w)] = l-T a T b uj 2 , (26) 
Im^(ju)} = (T a + T b )uj. (27) 

Note that 4>(ju>) (u> £ E) draws a parabolic curve on the complex plane. Proposition 
1 implies that w is the minimum positive solution of aTg((j)(jzu)) = ir/N — tut. Thus, 
Ke[(f>(jzu)} = Acos(ir/N) and lm[<p(jvj)] = Asin(7r/A^) for some A > 0. Substituting these 
into (|26|) and (|27l) . and eliminating A, we have the following equation. 

T a T b zu 2 + (T a + T b ) cot (j- - wrj vj - 1 = 0, (28) 

from which (|25p is immediately obtained. □ 

The oscillation frequency is analytically predicted in this theorem. Thus, it is possible to 
interpret the relation between the parameters and the frequency. Note that (|25|) is written 
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only in terms of the essential parameters shown in Table [T] We will demonstrate biological 
insights obtained from Theorem 1 in Section [6] 

Remark 2. The minimum positive solution of (|25|) can be efficiently obtained by bisection 
search for m G [0,7r/iW], because the right-hand side of ([25]) monotonically decreases for 
tu G [0, n/Nr], and the minimum positive solution of (|25|) exists in this region. It should be 
noted that w is obtained without the numerical computation as 




when the time delays are not considered, i.e., r ri — and r Pi = (i = 1, 2, • • • , N). 



4.2 The phase of oscillations 

The phase profile is determined from the eigenvector y* of UK\ in (|20b[) . In what follows, 
the phase profile is explored in an analytic approach. 

We first show that UK\ can always be transformed into a circulant matrix ( Davisl 19791) 
by similarity transformation. 

Lemma 3. Let D := diag(<ii, e?2, • • • , (In) G C NxN be defined by 



(-If- 1 i±*=l5*r (if n is odd) 



di~ < 



rifc=i a 



(30) 



-e J « 



(if AT is even) 



with 4* : = ^*(^*_i) V*-i)- Then, D 1 (UK)D = V, where V is the circulant matrix of the 
form 



V: 



n?=i£fcl*<yc(-l,-l,--- ,-1) (if ^ is odd) 
Ilfe=i ^fcl "^cyc(eTv , e^v , • • • , eiv ) (if JV is even). 



Since circulant ma trices are known to be diagonalized by the discrete Fourier transform 
matrix ( Davisl . 19791) . the eigenvector q := [qi, q 2 , ■ ■ ■ ,qN] T of V associated with the eigen- 
value <t)(jvj)e^ T is easily obtained as 

_j(-iye-^ (if N is odd) 

(if N is even). (6 ) 

Therefore, the phasor y — [yie^ 1 , y2&'' P2 , • ■ ■ , yNe : ' VN ] T is calculated from y = Dq. Finally, 
computing ifi = ipi + vjT pi (i = 1,2, •• • , N) yields the following analytic estimate of the 
phase. 

Theorem 2. Consider the gene regulatory networks modeled by fip. The phase shift (<Pi-\-i — 
ipi) between the (i + l)-th and the i-th protein is expected as 

<Pi+i -<P%= (Zi - — j 7r - wAt { . (32) 
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for i = 1, 2, • • • , N , where 



An := (r n+1 + r Pi ) - t, 



:= 



1 if 5, 



ifc5 i+1 



-1 
+ 1 



(33) 
(34) 



and w is given by Theorem 1. 

This theorem analytically predicts the phase difference between protein species. The 
constants Ari (i — 1, 2, • • • , N) are the discrepancy of the time delay of the i-th gene from 
the average delay r. The interpretation of Theorem 2 is given in Section [BJ 

Remark 3. When the delays are homogeneous between genes, i.e., r ri = r r2 = 



= T r 



and t, 



Pi 



I p 2 



'Pn ; 



32|) is independent of w, and it depends only on the activation- 



repression pattern of the regulatory network. 

4.3 The bias and amplitude of oscillations 

The bias and amplitude profiles can be predicted from x and \y\ in (j 15[) . 

We first show that the eigenvalues and eigenvectors of ICq(x, \y\) are obtained in a similar 
fashion to Lemma 1. 



Lemma 4 . For any given (x, y), the eigenvalues 
are given by 



N 



JJ T] k (x k -i,yk-i) 



fc=i 



1, 2, • • • , N) of the matrix ICo(x, \y\ 



(35) 



and the eigenvector associated with /xi is [771,771772, iii»=-i 

This lemma shows that the eigenvalues of Kq are located on a circle with center at the 
origin, but the angular position is different from that of UK\ (see Lemma 1). Following the 
similar argument as Lemma 1 and Proposition 1, we can show that [i\ should become the 
eigenvalue that corresponds to (f>(0) (— 1) in (|20ap . Therefore, the bias x* is obtained from 
the eigenvector of K$ associated with /ji, and we have the following theorem. 

Theorem 3. Consider the gene regulatory networks modeled by {IJj. The bias x and 
amplitude \y\ of oscillatory protein concentrations pt(t) (i = 1,2, ••• , N) are expected to 
satisfy both \36a\) and 



L36b\) simultaneously 

Xi+l 



21 i 

Vi+i 



Ci+i(xi,yj) 
>Vk-i)\ : 



^ i+1 (x lJ y l ) 



(36a) 
(36b) 



where w is given by Theorem 1. 

This theorem provides the equations that the bias and amplitude of oscillations should 
satisfy. Note that Xi and yi (i = 1, 2, • • • , N) are determined from the recursive equations 
(|36|) once (xi,yi) is determined. It also follows that 



N 



N 



Y[€i(xi-i,yi-i) 



(37a) 



(37b) 



11 



from P0|) . Thus, Xi and j/j (i = 1,2, ••• , N) could be obtained by numerically searching 
(2:1,2/1) so that (|37|) is satisfied under the constraint ([36]). However, the construction of a 
fast and reliable algorithm to compute x and \y\ remains a future challenge. 



5 Applications to Biochemical Oscillators 



In this section, we demonstrate the pro posed analysis me thod with two biochemical net- 
works. We first consider the Pcntilator (Tsai et al. . 120081) . a conceptual biochemical net- 
work consisting of N = 5 genes, to illustrate that the developed result is useful even for the 
large network. Then ; we investigate the oscillations of an existing oscillatory reactions, the 
somitogenesis clock (jHirata et all 120041 ) . 



5.1 Pentilator with time delay 

Pcntilator ( Tsai et al. , 2008) is a model of the gene regulatory network that is composed of 



N = 5 genes interacting in a cyclic way. Although the original model assumed that all the 
interactions are repressive, we here replace two of the repressors with activators in order to 
describe the effect of activators on the oscillation profiles. 

We consider the gene regulatory network depicted in Fig. @] (Top) . The dynamical model 
of the Pentilator in Fig. @] (Top) is then written for i = 1, 2, • • • , 5 as 



Pi(t) 



-an(t) 
cr t (t - 7 



i)-bpi(t), 



-J), 



(38) 



where /i(-) = / 2 (-) = /»(•) = F R (-) and / 3 (-) = / 5 (-) = F A (-). Note that time delays of 
tr anscription and tr anslation are explicitly modeled in (f38f , though they were not introduced 



Tsai et all ([20081 ) 



The parameters are set as follows: the degradation rates are a = 2.0[min~ 1 ] and b = 
0.2[min _1 ], and the synthesis rates are c = 0.3[min ] and /3 = 10[min ]. The Hill co- 
efficient is v = 2.0, and the time delays are t t := [1.8, 1.4, 1.1, 0.7, 1.0] T [min] and t p := 
[1.0, 0.8, 0.4, 0.4, 0.4] T [min], where the i-th entry of x r and t p is defined as T Ti and r Pi , 
respectively. 

From the above parameters, we have 



Q = 0.575, t = 1.8, T A = 1.1, 
At = [1.0, 0.4, -0.3, -0.7, -0.4] T , 



(39) 



where At is the vector whose z-th entry is At^ defined in (|33[) . Then, 4>{s)e ST in (|20[) becomes 
4>{s)e ST = (5s + l)(0.5s-|-l)e 1 - 8s , and the problem reduces to finding (ru,x,y) satisfying the 
conditions (CI) and (C2). 

The existence o f oscil lations for this system can be confirmed from the conditions presented 
111 iTakada et all (|2010l ). thus we hereafter analyze the waveform of the oscillations. Using 
Theorem 1, we first obtain the oscillation frequency as w = 8.98 x 10~ 2 [rad/min]. The 
phase profiles can be computed from Theorem 2, where Z\ = —l,Z-2 = —1,Z% = +1, Z4 = 
— l,Z$ = +1. The results are summarized in Table [2] Note that we define ipi = without 
loss of generality. 

To confirm these results, we conducted a numerical simulation of (|38[) . The time course 
of protein concentrations is presented in Fig. 0] (Bottom). The frequency and phase of 
the simulated oscillations are shown in Table [5] We see that the proposed analytic results 
provide approximate oscillation profiles with high accuracy. 
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Figure 4: (Top) Schematic diagram of the rearranged Pentilator considered in this example. 
(Bottom) Time course of the oscillatory protein concentrations in the rearranged Pentilator. 



Table 2: The estimated/actual frequency and phase of the oscillations 



Theorem 1 


8.98 xlCT 2 [rad/min] 


Simulation 


8.61 xlCT 2 [rad/min] 


Relative error 


4.30 % 





P2 


P3 


P4 


P5 


Theorem 2 [deg] 
Simulation [deg] 


146.1 
141.1 


108.5 
110.1 


255.6 
251.7 


218.1 
219.8 



5.2 Somitogenesis clock 

Somitogenesis is the process by which the somites of living organisms are produced. Recent 
biological studies revealed that somite segmentation of mouse embryos occurs every two 
hours, and the segmentation clock is primarily controlled by oscillatory gene expression of 
Hes7 gene (see Hirata et al. ( 20041) and references therein). 

The Hes 7 transcription is controlled by the self-negative feedback of Hes7 protein. Thus, 
the dynamics of Hes7 protein concentration can be described as 



h(t) = -an(t) + 



(40) 



Pi(t) = -bpt{t) + cri(t - r ri ). 



Note that this model is a special case of ([]} with N = 1. 

Theoretical studies of the self-negative feedback system ([¥0]) pre dicted that sh ort half- 
life of Hes7 i s a ke y to generate the oscillations of two-hour period (jMonkl . 120031) . Later, 
Hirata et al. (2004) conducted an experiment using Hes7 mutants with long half-lives and 
normal repressor activity, and confirmed that the mutants with long half-lives do not exhibit 
oscillations. Although the oscillations' existen ce was focused in many existing studies in- 



cluding the authors' work ( Takada et al. . 2010l ). the underlying mechanism that determines 



the profiles of oscillations are still unclear. 

In what follows, we demonstrate that the frequency profile can be analyzed with Theorem 
1 , then provide insig hts on the relation between the parameters and the frequency of the os- 
cillations. Following iHirata et al] (|2004f ). we here define the mRNA and protein degradation 
rates a and b as 



In 2 



In 2 



(41) 



where t r and t p denote mRNA and protein half-life time. The parameter values of wild-type 
Hes7 is also taken from lHirata et al.l ( 2004T ): t p = 20 [min],£ r = 3 [min],a = 0.231 [min _1 ],6 = 
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Figure 5: Time trajectory of Hes7 protein concentrations. Hes7 exhibits periodic oscillations 
with the two-hour period. 

0.0347 [min -1 ], c = 4.5 [min _1 ],/3 = 0.825 [min -1 ], r p = 17 [min],r r = 20 [min], 2/ = 2. Note 
that mRNA and protein concentrations are nor malized by the half maximal effective con- 



centration of the protein, or the constant p cr it in iHirata et al.l (|2004l) . 

Theorem 1 implies that the frequency profile is essentially determined from Q,t and 
Ta, whose meanings are listed in Table [TJ Thus, these quantities are calculated from the 
parameters as 

Q = 0.674, t = 2.23, T A = 16.6. (42) 

We obtain the frequency vu from Theorem 1 as w — 0.523 [rad/min], from which the period 
is calculated as 

— = 120.1 [min]. (43) 
w 

The numerical simulation of Hes7 concentrations indeed exhibits two-hour oscillations as 
shown in Fig. [5j We see that Theorem 1 successfully predicted the frequency of oscillations 
without simulating (|4"0"1) . 

Remark 4. Hirata et al.l (|2004f ) showed that the period of oscillations is approximately 



given by 2(r r +t p + ln(2)/a + ln(2)/6) with a and b defined by (|4ip . This estimate provides 
12 0.0 [min] in the abov e example. A more detailed comparison of Theorem 1 and the estimate 
in IHirata et al. ( 20041 ) is shown in Appendix [EJ It should be emphasized that Theorem 1 is 



useful for negative cyclic feedback networks consisting of any number of genes, while Hirata's 
estimate is applicable only for the self-repression case, i.e., N = 1. 

Theorem 1 also allows us to obtain more general insights on the frequency profile. The 
equation (|25[) implies that the parameters r and Q are important for characterizing the 
frequency. Based on Theorem 1, the relation between these parameters and the period 
2ir/m of oscillations is illustrated in Fig. [SJ Note that Q stands for the discrepancy of 
mRNA and protein degradation rates, and Q = 1.0 implies a = b = 0.0604[min~ 1 ]. We see 
that the frequency is sensitive to the sum of transcription and translation delay r — r r + t p 
rather than Q. Thus, the transcription and translation delays play an important role in 
regulating the two-hour cycle of the segmentation clock. 



6 Biological Insights 

In this section, we provide general biological insight based on Theorem 1 and Theorem 2. 
Specifically, we illustrate how the profiles of oscillations depend on parameters and structures 
of the network. 



14 




Transcription and translation time delay x [min] 

Figure 6: The period of Hes7 oscillations predicted from Theorem 1. The mean degradation 
rate of mRNA and proteins is set as Ta = 16.6, and R = 21.5, which is equivalent to the 
parameters in The oscillations illustrated in Fig. [5] correspond to the red point with 

the caption "Hes7 oscillation example." 

Frequency: The frequency profile is given by ([25j) in Theorem 1. Let w be defined by 
vj := vjTa with Ta in (j8|). From the definition, vo is a frequency normalized by the mean 
degradation rates of mRNA and proteins. Then, we can rewrite (|25|) as 

^=^2 (^cot 2 (|r^)+3 2 - cot (|r^)) - (44) 

where f is the time delay normalized by Ta (see © for the definition). This implies that time 
can be normalized by the mean degradation rates Ta without loss of generality. Therefore, 
the frequency profile can essentially be captured by the dimensionless quantities that appear 
in (|44]). 

We see that w depends only on the normalized time delay f , the number of genes TV and 
the discrepancy of the degradation rates Q (see Table [1] for the biological meanings). It is 
worth noting that as a result of cyclic feedback, the mean time delay over all genes, f, plays 
a decisive role rather than the individual delays at each gene. Thus, the frequency can be 
less sensitive to the variation of the delays in intermediate reactions. 

Looking more closely into (|44j) . we obtain more quantitative relations between these pa- 
rameters and the frequency as follows: the normalized frequency w monotonically increases 
as 

(i) the average time delay of transcription and translation process (f ) decreases 

(ii) the number of genes (N) decreases 

(iii) the mRNA and protein degradation time gets close to each other. 

These insights are confirmed by numerical simulations as shown in Fig. [7] Note that the 
vertical axis is the period of oscillations, where the time is normalized by Ta- We see that 
the solution of (JUJ) successfully predicts the oscillation profiles. 
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Figure 7: The normalized period 2tt/vd of oscillations in terms of f. 
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Figure 8: The phase profiles of oscillations for various f. When At 
intact, while the frequency changes drastically. 



0, the phase stays 



In particular, it is observed that the period is almost linearl y proportiona l to t he mean de- 
lay f. This observation is essentially consistent with the ones in Monk ( 2003 ) and Wang et al.l 



(|2004f ) , where the period of oscillations was shown to be proportional to the total time delay 



^^ =1 (r ri + r Pi ) based on numerical simulations. In fact, we have 



-^47v+25> ri +T p J — . (45) 

by approximating (|4"l)l with \J\ + x 2 ~ l+ir 2 /2 and cot(x) ~ x" 1 , and the normalized period 
of oscillations 2tt/w is indeed proportional to the sum of delays in Fig. [7] (see Appendix [C] 
for the details of the approximation). 

Remark 5. The frequency profile obtained from (pFfj) becomes less accurate, when the 
parameters R and v are so large that the oscillation waveforms are saturated, because the 
analysis is based on ([T0| . However, the authors have observed by numerical simulations that 
the qualitative insights obtained above hold, even when the waveform is much distorted from 

dia. 



Phase: The phase of oscillations is given by (|32|) . Since zuAt = tzjAt with Af, := At;/Ta, 
(|3"2"j) can be rewritten with the dimensionless quantities shown in Table Q] We see that the 
phase is determined from the activation and repression pattern Zi of the network and otA-t^. 



1G 



The constants Afj are the discrepancy of the delay from the mean delay r. Thus, the 
delays have no impact on the phase profile if they are homogeneous, while they do have an 
impact for the frequency profile. An example is illustrated in Fig. where At = 0. We see 
that only the frequency is affected by the difference of the mean delay f . 

The phase profile can be summarized as follows. 

(i) The difference of the individual time delays from the average f affects the phase of 
oscillations. 

(ii) The repression/activation of transcription causes phase delay/lead of the following 
protein oscillations. 

7 Conclusion 

We have developed a systematic method to explore the profiles of oscillatory protein concen- 
trations in biochemical networks with negative cyclic feedback. First, we have formulated 
harmonic balance for biochemical oscillators. Then, the relation between the reaction rates 
and the frequency, phase and amplitude of the oscillations has been analytically obtained. 
The proposed method has been demonstrated with the Pentilator and the somitogenesis 
clock, and the analysis result has shown quantitative agreement with the existing biological 
experiments. 

A distinctive feature of the proposed approach is that the oscillation profiles are written 
in an analytic way. Thus, we can obtain detailed biological insight, which provides use- 
ful guidance in synthetic gene circuit design. We have revealed dimensionless parameters 
that primarily encode the profiles of oscillations, and analyzed the qualitative properties of 
frequency and phase in terms of these parameters. 
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A The proof of Lemma 2 

The poles of %\{s) are obtained by solving 



\<j>{s)e ST I - UK X \ = \<f>{s)e ST I - A| = 0, 
where A := diag(Ai, A 2 , • • • , X N ) € C NxN . It follows that 



(46) 



N 



\<t>{s)e ST I - A| = l[Ws)e ST - Ai)(=: *(«)). 



(47) 
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Figure 9: Comparison of the period of oscillations estimated from Theorem 1. iHirata et al 
(|2004f ) and the simulation of (|40|) . 



DefineX(s) :— J[ i=1 i=jLl e ,((/)(s)e ST — A*). Then, we can write x(s) = (</>(s)e sr — A£)(<75(s)e ST — 
Xt')X(s). It should be noted that Hi(s) is a retarded time delay system, thus the dominant 
pole is located at the rightmost position in the complex plane. 

In what follows, we prove the sufficient condition. It follows from the definition of f2+ 
that the roots of X(s) = lie in the open left half complex plane. On the other hand, 
4>(s)e ST — X( = and (j>(s)e ST — Ag> = have at least one root at s = jujQ and s = —juio, 
respectively, and the rest in the open left half plane. In particular, the roots s — ±jcoo of 
these equations are single roots, because uj satisfies \cj){ju)o)eP' JJ ° T \ = \4>(joJo)\ — \Xe\ and 
|</3(jo;)| monotonically increases with respect to lj > 0. Therefore, 'Hi(s) has a pair of poles 
at s = ±jwo with multiplicity one, and the rest in the open left half plane. 

The proof of the necessary condition is omitted, since it can be easily obtained from the 
above discussion. □ 



B 



Comp arison of Theorem 1 and the estimate in lHirata et al. 

hood ) 



In IHirata et al. the period of oscillations is analytically estimated as 2(r r + t p + 

ln(2)/a + ln(2)/6 . where a and b are defined by (HT]) . We here compare the period obtained 
from Theorem 1. IHirata et al.l ( 2004 ) and numerical simulation of (|40| . The result is shown 
in Fig. 1] 

We see that Theorem 1 and IHirata et al. ( 2004 ) provide close estimation for a wide range 
of time delay. The actual oscillation period starts to deviate from these estimations as the 
time delay increases, because the waveform tends to be distorted. Nevertheless, the relative 
error at r = 50 [min] is —6.28% for Theorem 1 and —9.18% for the Hirata's estimate. Thus, 
these estimates can be useful for gaining qualitative insight of the oscillation period. Note 
that Theorem 1 is us eful for bioch e mical networks consisting of a large number of genes, 
while the estimate in IHirata et al.l (|2004J ) is applicable only to the self-negative feedback 
case, i.e. N = 1. 
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C Approximation of (1441) 



The right-hand side of (|44|) is rewritten as 

± {cot (£ - a? ) - 1) } . 

As t -> oo, ir/N — wt — > follows, thus Q 2 tan 2 (7r/iV — rof ) becomes sufficiently small so 
that we can approximate the square root by y/l + x 2 ~ f + x 2 /2. Then, we have 

w - tan - wr^j . (48) 

When ir/N — wf is close to zero, tan(7r/7V — wf) ~ ir/N — wf. Thus, w ~ l/2(n/N — wt), 
and the normalized period of oscillations is obtained as 

2 W j 

— ~4iV + 2iVf = 4JV + 2V(r ri +r Pi ) — . (49) 
w r—t J-A 
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